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COMPARISON OF RESPONSE SURFACE AND KRIGING MODELS IN THE 
MULTIDISCIPLINARY DESIGN OF AN AEROSPIKE NOZZLE 

TIMOTHY W. SIMPSON' 

Abstract. Tlie use of response surface models and kriging models are compared for approximating 
non-random, deterministic computer analyses. After discussing the traditional res}x)iise surfa(*e approach 
for constructing polynomial models for approximation, kriging is presented as an alternative statistical- 
based ai>proximation method for the design and analysis of computer experiments. Both approximation 
metho<ls are applied to the multidisciplinary design and analysis of an aerospike nozzle which consists of a 
computational fluid dynamics mcxlel and a finite element analysis nuxlel Error analysis of the response 
surface and kriging models is performed along with a graphical comparison of Uie approximations. Four 
optimization problenns are formulated and solved using lx)th approximation models. While neither 
approximation technique consistently outperforms the other in this example, the kriging models — using 
only a constant for the underlying global model and a (jaussian correlation function- jx'r for in as well as 
the second order ix>lynomial response surface modt^ls. 

Key words. response surface nuxlels, kriging, multidisciplinary design 

Subject classification. Applied and Numeri(*al Methods 

1. Introduction. Current engineering analyses rely heavily on running complex, and often 
exj>ensive, computer analysis codes. Despite the steady and continuing growTh of computing power and 
speed, the complexity of these codes seems to keep pace with computing advances. Statistical techniques 
are wudely used in engineering design to construct approxhuHtions of these analysis codes; these 
approximations are then used in lieu of the actual analysis codes, offering the following l^enefits. 

• They yield insight into the relationship IxtW'een (output) responses, y, and (input) design variables, 

X. 

• They provide fast analysis tools for optimization and design space exploration since the cheaj>to-run 
approximations replace t he expeiisive-to-run computer analyses. 

• They facilitate the integration of discipline dependent analysis codes. 

The most common method for building approximations of computer analyses is to apply design of 
experiments (DOE), response surface (RS) models, and regression analysis to build second order 
jx)lynomial approximations of the computationally expensive analyses. For example, the Robust C^oncept 
Exploration Method (see, e.g., [5] and l6]) has been developed to facilitate quick evaluation of different 
design alternatives, identify important ck^sign drivers, and generat<^ robust, top-level design spxicifications 
using DOE, RS nuxiels, and the compromise Decision Support Problem [25]; it has been successfully 
applied to the multiobjective design of a high speed civil transport (s(^e, e.g., [5] and [19]). a family of 
General Aviation aircraft [40], a turbine lift engine [I8j, and a flywheel [22]. In other work, the Variable 
Complexity Response Surface Mcxleling (VCRSM) methcxl (see, e.g., [12] and [13]) uses analyses of 
varying fidelity to reduce the design space to the region of interest and build response surface models of 
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incretisiiig accuracy. The VCRSM method employs DOE and RS modeling techniques and has been 
successfully applied to the multidisciplinary wing design of a high speed civil transport (see, e.g., [14] and 
[17]). to the analysis and design of composite curv^ed channel frames [24], and to reduce numerical noise 
inherent in structural analyses (see, e.g., [14] and [45]) and shape design problems using fluid flow analysis 
[29;. A review of several applications of response surfaces in mechanical and aerospace engineering design 
is given in [41;. 

The use of response surfaces for approximating deterministic (i.e., non-random) computer analyses is 
statistically questionable due to the lack of random error in the computer model (cf., [41J). A more 
appropriate and perhaps more “statistically sound” method for approximating deterministic computer 
experiments is through the use of kriging [9] models which are also referred to as the Design and Analysis 
of Computer Exix^riments (DACE) models (see, e.g., [4], [20], and [38]). The validity of the kriging model 
is not dependent on the existence of random error and may therefore be l^etter suited for applications 
involving determiiiLstic computer experiments. Furthermore, kriging models interpolate between data 
ix)ints which may le yield more accurate results since computer experiments typically do not contain 
random error (i.<‘.. you get the same output when you use the same input). 

Booker 2 contrasts traditional DOE and RS modeling with DACE models. In the “classical” design 
and analysis of j)hysical ('Xieriments, random variation is accounted for by spreading the sample points 
out in the d<*sign spac<‘ and by taking multiple data points (replicates), see FiC. 1. Sacks, et al. [38] state 
that the LichLssk*al notions of experimental blocking, replication, and randomization are irrelevant when 
it comes to determiiiLstic computer experiments; thus, sample points should l)e chosen to fill the design 
space. They sugg<>st minimizing the integrated mean squared error (IMSE) over the design region by 
using IMSE-opt imal dc'sigiis such as the one shown in the top right corner of FiC. 1. 
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As shown in the middle of FiG. 1, response surface modeling typically employs least squares 
regression to fit a polynomial model to the sampled data; kriging models are chosen to interix)late the 
data and are fit using maximum likelihood estimation (see. e.g., [20]). Validation of RS models is based 
on: (a) testing statistical hypothesis (t-tests and F-statistics) derived from error estimates of the 
variability in the data, (b) plotting and checking the residuals, and (c) computing the ratio of the 
nuxlel sum of squares to the total sum of squares (see, e.g., [28]). Sacks, e( al. [38] and Welch, et al. [47] 
both state that statistical testing is inappropriate when it comes to deterministic computer exj>erimenls 
which lack random error; cross-validation and integrate mean square error can be utilized to assess the 
accuracy of a kriging model. 

DACE and kriging models have found limited application in engineering design perhai)s because of 
the lack of readily available software to fit kriging models, the added complexity of fitting a kriging 
model, or the additional effort required to use a kriging model. To clarify this last p>oint, prediction with 
a kriging model requires the inversion and multiplication of several matrices, and the kriging mexUd ckx^s 
not exist as a ‘‘closed- for in'* ix)lynomial equation, RS nKxlel prediction only requires computation of 
simple j>olynomial equation onc(^ a model has been fit. The goal in this pajx^r is to examine the added 
computational expense required to perform kriging and compare the predictive capability of kriging and 
RS models. 

In Section 2 an overview of the statistical and mathematical foundations of response surface modeling 
and kriging is given. In Section 3 the multidisciplinary design of an aerospike nozzle is introduced; it 
serves as a test problem to compare RS and kriging models for approximation. In Section 4 the RS and 
kriging models are constructed for the aerospike nozzle example. In Section 5 four optimization problems 
are formulated and solved using the RS and kriging models; Section 0 contains a discussion of ongoing 
work. 

2. Statistical Approximations for Computer Experiments. Building approximations of 
computer analyses typically involves: (a) choosing an experimental design to sample the computer 
analysis code, (b) choosing a nuxlel to represent the data, and (c) fitting the model to the observed data. 
There are a variety of oi)tioiLs for each of these steps as shown in FiG. 2, and some of the more prevalent 
ones have been highlighted. 


3 




Fig. 2. Techniques forbuildingapproximations\4l\. 


Simpson, et al. [41j discuss several of the advantages and disadvantages of the highlighted 
approaches listed in FiG. 2, namely, resp>onse surface methodology, neural networks, inductive learning 
and kriging. In this work, the model choice and model fitting portions of building approximations are the 
primary concern; responst^ surface models are discussed in Section 2.1 and kriging models in Section 2.2. 
It is assumed that the reader has some knowledge of DOE and RS modeling and little to no knowledge of 
kriging. 

2.1. Overview of Ftesponse Surface IVlodeling. Response surface modeling techniques were 
originally develoj>ed to analyze the results of physical experiments and create emi)irically-based models of 
the observed resjxjnsc^ values. Resix)nse surface modeling jx>stulates a model of the form: 

(1) y(x) - f(x) ^ £ 

where y(x) is the unknown function of interest, f(x) is a known polynomial function of x, and e is random 
error which is assumed to \ye normally distributed with mean zero and variance The individual errors, 
£j, at each ol)servation are also assumed to l)e indejxmdent and identically distributed, llie ix)lynoniial 
function. f(x), used to approximate y(x) is typically a low order polynomial which is assumed to Ix^ either 
linear iis given by Eqn. (2): 


(2) y=3o + XPiXi 

i^l 

or (piadratic as given by Eqn. (3): 


(3) 


k k 

y=Po + ZPiXi + XPixf + 

i = ! 


^ ^ Pi] 

i<j 
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The parameters pQ, pj, p^^, and p^j, of the polynoiniaLs in Eqns. (2) and (3) are determined through 
least squares regression which minimizes the sum of the squares of the deviations of the predicted values, 
y (x), from the actual values, y(x). In order to fit the model to the observed data using least squares 
regression, the coefficients of Eqns. (2) and (3) can be estimat'd using IOQN. (4). 

(4) p-[X’Xi-iX’y 

In Eqn. (4), X is the design matrix of sampled }X)ints, and y is a column vector containing the 
corresjx)iiding values of the re^sponse. Tor more details on least squares regression or polynomial RS 
modeling see, e.g., [28]. 

2.2. Overview of Kriging. 

2.2.1. Mathematics of Kriging, Kriging, or DA(3^ modeling, postulates a combination of a 
polynomial model plus dei)artures of tlu^ form given by Eqn. (5): 

(5) y(x) = f(x) - Z(x) 

where y(x) is the unknown function of interest, f(x) Is a known polynomial function of x, and Z(x) is the 
realization of a normally distributed CJaussian random process with mean zero, variance o^, and non-zero 
covariance. The f(x) term in Eqn. (5) is similar to the jx>lynomial model in a response surface and 
provides a ‘'globaf' nuxlel of the ck^sign space; in many cases f(x) is simply taken to be a constant term p 
(see, e.g., [20], [38], and [47]). 

While f(x) *‘gloljally" approxiniatc^s the design space, Z(x) creates ‘dcx'alized" deviations so that the 
kriging model inter|K>lates the n^, sampled data points. The covariance matrix of Z(x) is given by Eqn. 

( 6 ) . 

(6) Cov[Z(x‘),Z(xJ)] - R([R(x^xJ)]) 

In Eqn. (6), R is the correlation matrix, and R(x*,xj) is the spatial correlation function IxUw^een any two 
of the iig sample points x* and xk R is a n,^ x iig symmetric, jx)sitive definite matrix with ones along the 
diagonal. The correlation function R(x',xj) is specified by the user; Sacks, et al. [38; and Kcxdiler and 
Owen [20] discuss several correlation functions. The Caussian correlation function in Eqn. (7) is 
employed in this w^ork. 

(7) R(x',xj) = exp(-Xk=i®k|xk -X'k| ) 

The 0]^ in EQN, (7) are the unknown parameters used to fit the model, n^jv is the number of design 
variables, and x^* and X|^j are the k*^ components of sample points x* and xh In some cases, using a 
single correlation parameter gives sufficiently good results (see, e.g., [4], [30], and [38]). 

Predicted estimates, y(x), of the response y(x) at untried values of x are given by Eqn. (8): 

(8) 9=p+rT(x)R-i(y-Fp) 
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where y is the column vector of length iig which contains the values of the response at each sample point, 
and F is a column vector of length n^ which is filled wdth ones when f(x) is taken as a constant. In EQN. 

(8) , r(x) is the correlation vector of length iig l^etween an untried x and the sampled data points [x^ x*^, 
x"«J and is given by Eqn. (9). 

(9) = [R(x,x^), R(x,x^), R(x,x“^)j"^ 

In Eqn. (8), P is estimated using Eqn. (10). 

(10) p =(FTR-iF)F^^R-iy 

The est imate of the variance, o , of the sample points from the underlying global model (not t he variance 
in the observed data itself) is given by EQN. (11): 



where f(x) is assumed to be the constant p . The maximum likelihood estimates (i.e., “best guesses") for 
the 01^ in Eqn. (7) used to fit the model are found by maximizing EqN. (12) over 0k > 0 (see, e.g., [4]). 


( 12 ) 


max 

e, >0 



Both and detR are functions of 0k- While any values for 0k create an interpolative model, the ”l3est" 
kriging ukkIcI is found by solving the k-diniensional unconstrained, non-linear, optimization problem given 
by Eqn. (12), 

2.2.2. Applications of DACE and Kriging. DACE and kriging models have found limited use in 
engineering design applications since its introduction into the literature by Sacks, et al. [38], Ciunta [11] 
has i^erformed a preliminary investigation into the use of DACE modeling for the multidisciplinary design 
optimization of a High Speed Civil Transport aircraft. He explores a 5 and a 10 variable design problem, 
olx>er\dng that the DACE and response surface modeling approaches yield similar results due to the 
(quadratic t rend of the res{X)iises. Osio and Amon [30] have developed an extension of DACE modeling for 
numerical opt imization which uses a multistage strategy for refining the accuracy of the model; they have 
applied their approach to the thermal design of an einlx^dded electronic package wdiich has 5 design 
variables. Welch, et al. [46] descrilx^ a kriging-based approximation methodology which they use to 
identify imixirtant variables, detect curvature and interactions, and produce a useful approximation model 
for two 20 variable problems using only 30-50 runs of the computer code; they claim their method can 
cojxi with up to 30-40 variables provided factor sparsity can l>e exploited. Booker, et al. [3] solve a 31 
variable helicopter rotor structural design problem using an approximation method based on kriging, 
Booker [2] extends the helicopter rotor design problem to include 56 structural variables to examine the 
aeroelastic and dynamic response of the rotor. Trosset and Torezon [44] have develojxid a numerical 
optimization strategy which incorporates DACE modeling and pattern search methods for global 
optimization. Cox and John [7] have developed the Sequential Design for Optimization method which 
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uses lower confidence bouncis on predicted values of the response for the sequential selection of evaluation 
points during optimization. Both approaches have shown improvements over traditional optimization 
approaches when applied to a variety of standard mathematical test problems. 

2,3. One Variable Example of Response Surface and Kriging Models. A simple one variable 
example^ bests illustrates the difference between the approximation ca|>abilities of a second order RS 
model and a kriging model. Su and Renaud [42] formulated this example to demonstrate some of the 
limitations of using second order RS models; see FiG. 3. They fit a second order response surface using 
least squares regression to five sanqile j>oints from a fabricated eighth order function within the region of 
the optimum (x = 932). A kriging model using a constant for the global model and the Gaussian 
correlation function of EqN. (7) is fit to the same five points; the original function, the five sample points, 
and the RS and kriging models are shown in Fl(i. 3. 

Immediately evident from 3 is fact that the kriging nuxlel interpolates the data ix>ints, 

approximating the original function l^etter than the response surface model and predicting an optimum 
wdiich is much closer to the actual optimum. It is imjx>rtant to notice that outside of the design space 
defined by the sample points (920 < x < 945), neither model predicts well as expected; the kriging model 
returns to the underlying global model which is a constant. This is typical Ixdmvior for a kriging nuxlel; 
far from the sample points, the kriging model returns to the underlying global model since t he influence' of 
the sample points has exix>nentially decayed away outside of the design space. 



Sixteen evenly spaced jx)ints (not including the sample points) are taken from within the range (920 
< x < 915) to check the accuracy of the two approximation models. The maximum alisolute error 
(difference between the actual and predicted values), the average absolute error, and the root mean square 
error (MSE), which is: 


(T3) 


root MSE = 


. _ Itt (Vi ■'f f 

"I %rror 
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where iierror is the number of points (= 16), are listed in TABLE 1. Based on the error analysis in TABLE 
1, it can be concluded that the kriging model approximates the original function Ijetter since it has a 
lower root MSE, average error, and maximum error* A more involved multidisciplinary design example is 
described in the next section. 


Table 1. Erroi'analysisforonfvariableex ample. 


Max A BS (error) 
Avg A BS (error) 
root MSE 


2nd Order RS Model 

Kriging Model 

37m 

2307 

1.911 

0.77(> 

2.155 

1.004 


3. Aerospike Nozzle Design Example. The multidisciplinary design of an aerospike nozzle has 
l>een selectetl as t he test problem for comparing the predictive capability of RS and kriging models. The 
linear aerospike rcx'ket engine is the propulsion system proix>sed for the VentureStar [4.3] reusable launch 
vehicle illustrated in FKL 4. 



KKj . 4. \ f-iit ui'f.Star7'eusabtfAauncttw:}tirleanthli‘ttf-:araei‘ospikfip7'op'ul»ionfiy:ttf^-in\2i\. 

Tlie aerospik<‘ rocket engine consists of a roc’kel thruster, cowl, aerospike nozzle, and plug base 
regions as shown in FlC. 5. The aerospike nozzle is a truncated spike or plug nozzle that adjusts to the 
ambient pressure and integrates well with launch vehicles [34]. The flow^ field structure changes 
dramatically from low altitude to high altitude on the spike surface and in the base flow region (cf., [15], 
[27], and [36]). Additional flow is injected in the base region to create an aerodynamic spike [16] which 
gives the cU'rospike nozzh' its name iind increases tin* base pressure and contribution of tli<^ base region to 
the aerospike thrust. 
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Fici. 5. Afro spikecompo7iP.rit sand flow fmldrh ararf prisfic s\2 \ \ , 


The analysis of the nozzle involves two disciplines — ^aerody nan lies and structures — since there is an 
interaction between the structural displacements of the nozzle surface and the pressures caused l)y the 
varying aerodynamic effects. Thrust and nozzle w'all pressure calculations are made using comj)utat ional 
fluid dynamics (CFD) analysis and are linked to a structural finite element analysis model for determining 
nozzle weight and structural integrity. A mission average engine sjx^cific impulse and engine 

thrust/weight ratio are calculated and used to estimat e vehicle gross-lift -off- weight (CLOW) based on 
data supplied by Rxx'ketdyne. The multidisciplinary domain decomposition is illustrated in FiG. (i. 
Korte. et al. [21] provide additional details on the aerodynamic and structural analyses for the a<^rospike 
nozzle. 



Fk J . 6 . Xfulfidiscipiiriary domain df^rmnposition [2 1 1 . 

For this study, three design variables are considered for the multidisciplinary design of the aerospike 
nozzle: thruster angle, base height, and length as shown in FiG. 7. The thruster angle is the entrance 
angle of the gas from the combustion chaml)er onto th(' nozzle surface; the base height and length refer to 
the solid portion of the nozzle itself. A quadratic model is created to generate values of spline knot 
surface angle slope and exit angle which define the nozzle profile, corresi)onding to sjx^cific values of 
thruster angle, base height, and length. 



Fig. 7. Nozzlegeorrietrydesig7iva7'iabl^:s\2l\. 
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Bounds for the design variables are set based on information from Boeing Rocketdyne and viable 
nozzle profiles from the quadratic model based on all combinations of thruster angle, height, and length 
within the design space. Second order RS models and kriging models are developed for each 
response — thrust, weight, and GLOW — in the next section; optimization of the aerospike nozzle using the 
RS and kriging models for different objective functions is performed in Section 5. 

4. Approximations for the Aerospike Nozzle Problem. The data used to fit the RS and 

kriging nuKlels is obtained from a 25 point orthogonal array (lx)se 5 3 I oarand) from [33j. The sample 
ix)ints are illustrated in FiG. 8 and are scaled to fit the three dimensional design sp>ace defined by the 
lx)unds on tlie thruster angle, base height, and length. 


Length 



Section 4,1 contains the resp>onse surface models which are fit to the data; Section 4.2 contains the 
kriging ukkIcIs. Error analysis of the response surface and kriging models is discussed in Section 4,3, and 
a graphical comparison of the approximation models is given in Section 4.4. 

4.1. Response Surface Models. The RS models for weight, thrust, and GLOW are obtained using 
ordinary least squares regression techniques and JMP [39]. The corresponding RS models are given in 
Eqns. (14)-(1()). The equations have l^een scaled against the baseline design due to the proprietary 
iiatun^ of some of the data. 

Weight = 0.810 - 0.116*a + 0.121*h + 0.152*1 -h 0.0()5*a^ - 0.025*a*h + 

(14) 0.0013*h^ - 0.0539*a*l - 0.0131*h*l + 0.0301*E 

Thrust = 0.9968 + 0.00031*a + 0.0019*h + 0.0060*1 - 0.00175*a^ + 0.00125*a*h - 

(15) 0.001 l*h^ + 0.00125*a*l - 0.00198*h*l - 0.00165*E 

GLOW = 0.9930 - 0.0270*a ^ 0.0065*h - 0.0265*1 + 0.0307*a2 - 0.0163*a*h + 

(16) 0.0100*h^ - 0.0226*a*l + 0.0151*h*l -r 0.0195*P 


The R^. R*^-adj listed, and root MSE values for each of these second order RS models are summarized 
in Tablk 2. As evidenced by the high R'^ and R/^-adjusted and low root MSE values, the second order 
ixjlynoniial mcKlel appears to capture a large portion of the observed variance. 


10 



T A B L E 2. M odeldiagnos ticsofrespon su7'facemodeLs\ 


Respoiise 


Measure 

Weight 

Thrust 

GLOW 

TP 

(TMT 


o:m 

R/^-adjustecl 

0.977 

0.996 

0.953 

root MSE 

1.12% 

0.01%, 

0.25% 


4.2. Kriging Models for the Aerospike Nozzle Problem. The krigiiig models are from 
the same 25 sample points used to fit the response surface models in Section 4,1. A constant term (i.e., 
the mean of the data) is selected for the underlying gloljal model, and the Gaussian correlation functiom 
EQN. (7), is utilized for the local departures determined by th<' correlation matrix R. 

Initial investigations revealed that a single 0 parameter Wris insufficient to accurately model the data 
due to sc*aling of the design variables. Therefore, a simple 3-D exhaustive grid search with a refinable 
step size was used to find the maximum likelihood estimates for the three 0 parameters needed to obtain 
the kriging model. The resulting maximum likelihood estimates for the three 0 parameters for the 

w'eight, thrust, and GLOW models are summarized in TABLE 3; these values are for the scaled sample 
ix)ints. 


Table 3. Tltetapararnetf^Tsforki'ig-ingmodeLs. 


Response 

Weight 

Thrust 

GLOW 

D34BTj 

(mr 

062 


0.50 

2.437 

0wh = 

0.65 

0.537 


With these parameters and the corresponding 25 sampk‘ points, the kriging models are fully 
specified, A new point is predicted using these 0 values in combination with Eqns. (8)-(10). The 
accuracy of the RS and kriging models is examined in the next two sections. 

4.3. Error Analysis of Response Surface and Kriging Models. An additional 25 randomly 
selected points are used to verify the accuracy of the RS and kriging models. Error is defined as the 
differenc’e between the actual response from the computer analysis, y(x), and the predicted vahie, y (x), 
from the RS or kriging model. The maximum al^sohite error, the average al^solute error, and the root 
MSE, Eqn. (13), for the 25 randomly selected points are summarized in TABLE 4, 


Table 4. Erroi'analysisofapproximationmodf^ls. 


Se condOrdevRe sponseSur f ace M odcls 


Weight 

Thrust 


Max ABS(error) 19.57% 

Avg ABS(error) 2.44%) 

root MSE 4.54% 

K rising M odels (wit ha 

o:o32%r 

0.012% 

0.015% 

mstantter'i 

08% 

0.53% 

0.90% 

’n) 


Thrust 

m 

Max A BS (error) 17.23% 

Avg A BS (error) 2.51% 

root MSE 4.37%, 

04148%: 

0.012% 

0.018% 


11 



For the weight and GLOW responses, the krigiiig models have lower maximum al>?olute errors and 
lower root MSEs than the RS models; how^ever, the average al^solute error is slightly larger for tlu^ kriging 
iiKKlels indicating that the average magnitude of the prediction error is larger for the kriging models than 
the RS models. As for thrust, the RS models are slightly better than the kriging models according to the 
values in the table; the maximum absolute error and root MSE are slightly less wdiile the average aljsolute 
errors are essentially the same. It is not surprising that the RS model predicts thrust better; it luis a very 
high value, 0.998, and low nx)t MSE, 0.01%. It is reassuring to note, however, that the kriging model, 
despite using a constant term for the underlying global model, is only slightly less accurate than the 
corresjx)iiding RS model. It appears that both approximations predict reasonably well wnth the kriging 
models luiving a slight overall advantage Ixicause of the lower root MSE values. 

4,4. Graphical Comparison of Response Surface and Kriging Models. In addition to the 
error analysis of Section 4.3. a graphical comparison of the RS and kriging models luis l^een i>erformed to 
visualize differences in the two approximation models. In FiGS. 9-11, 3-D contour plots of thrust, weight, 
and GLOW as a function of angle, length, and base height are given. In each figure, the same contour 
levels are used for th<‘ RS and kriging models so that the shapes of the contours can he compared. 



Kk^ . 9. Rf.sponHesurfacfiandki'igingrnodels forihrust . 

In Fig. 9, it can Ix' seen that the contours of thrust for the RS and kriging models are very similar. 
As evidenced by the high R^ and low root MSE values, the RS model fits the data quite well, and it is 
reassuring to note that the kriging model resembles the RS model even through the underlying global 
model for the kriging nKxlel Ls just a constant term. 
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Fig . 10. Responsesurfaceandkriffingmodelsfoi'U'eight. 

The contours of the RS and kriging models for weight in FlCi. 10 are also very similar, but the 
influence of tlu^ hK alized perturbations caused by the Gaussian correlation function in the kriging mo<lel 
can begin to Ix' s<k'ii. The error analysis from Section 4.3 indicated that the kriging model for weight is 
slightly mor(' ficcural<' than the RS model which may result from the small localized variations. 



A. Ismrif^irirview. B. Endview. 

P'lG. 1 1. Responsesurfaceandh'ig^mgTnod€:lsfoi'GLOW. 


The general shape of the GLOW contours is the same in FiG. 11 A, however, the size and sha|x^ of the 
different contours, particularly along the length axis, are quite different. The end view along the length 
axis in FiG. llB further highlights the differences Ix^tweeii the two models. Notice also in FiG. llB that 
the kriging nuxlel predicts a minimum GLOW within the design space centered around Height=-0.8, 
Angle— 0, along the axis defined by 0.2 < Length < 0.8; this minimum was verified through additional 
experiments. 

5. Optimization Using Response Surface and Kriging Models. The true test of the accuracy 
of the RS and kriging models comes when the approximations are used in optimization. It is crucial that 
any approximations used in optimization prove reasonably ac'curate, lest they lead the optimization 
algorithm into regions of bad designs. Trust Region approaches (see, e.g., [23] and [35]) and the Mcxlel 
Management framework (see, e.g., 4] and [10]) have l)een develojxxl to ensure that optimization 


13 



algorithms are not led astray by inaccurate approximations. In this work, however, the focus has Ix^en on 
developing the approximation models, particularly the kriging models, and not on the optimization itself. 

Four different optimization ])roblems are formulated and solved to compare the accuracy of the RS 
and kriging models: (1) maximize thrust, (2) minimize weight, (3) minimize CLOW, and (4) maximize 
thrust/weight ratio. The first two objective functions represent traditional single objective, single 
discipline optimization problems. The second tw'* objective functioiLs are more characteristic of 
multidiscipliiiary optimization; minimizing CLOW or maximizing the thrust/weight ratio requires trade- 
offs Ix'tween the at^rodynaniics and structures disciplines. For each objective function, constraint limits 
are placed on the remaining responses; for instance, constraints are placed on the maximum allowable 
weight aiul CLOW and the minimum allowable thrust/weight ratio when maximizing thrust. However, 
none of the constraints are active or binding in any of the final results. 

Each optimization problem is solved using: (a) the RS model approximations and (b) the kriging 
iiKxlel approximations for thrust, weight, and (JLOW. The optimization is ix^rformed using a Ceneralized 
Reduced Cradient (CRC) algorithm. Three different starting points are used for each objective function 
(the lower, middle, and upper bounds of the design variables) to assess the average number of analysis 
and gradient calls necessary to obtain the optimum design within the given design sj>ace; the solutions for 
each objf^ctive for each approximation converge to the same oi)timuni despite the initial starting point. 
The same paraiiK'ters (i.e., step size, tolerance, constraint violation, etc.) are used within the CRC 
algorithm for ecich optimization. The optimization results are summarized in TABLK 5. Design variable 
and resjx)ns<" values have Ix^en scaled as a percentage of tlu‘ l>aseliiie design due to the proprietary’ nat un^ 
of some of the data. 


Table 5. OptimizationresuU susingvespmisesur farf^^a^idkriginginodels. 



Avg. # of 
Analysis 
Calls 

Avg, # of 
Cradient 
Calls 

Optimum Design 

Predicted Optimum 

Verified 

Optimum 

% Error* 

Maximize Thrust 




Angle 


Thrust 

1.0016 

1.0013 


RS 

27 

4 

Height 

-0.433 

Weight 

0.9450 

0.9476 

-0.27% 

ModeAs 



l.eugtb 

1.000 

Tln-/Wt 

1.0141 

1.0134 

0.07% 






GLOW 

0.9724 

0.9759 

-0..16% 




Angle 

0.656 

4'hrust 

1.0016 

1.0014 

0.02% 

KiAffing 

62 

,5 

Height 

-0.627 

Weight 

0.9385 

0.9155 

2.51% 

ModeAti 



Length 

1.000 

Thr/Wt 

1.0157 

1.0210 

-0.51% 


j 

1 

i 




G1.0W 

0.0690 

0.9683 

0.08% 


Minimize Weight 





Angle 


Thrust 

0.9957 

0.9957 


RS 

29 

3 

Height 

-1.000 

Weight 

0.7584 

0.7496 

1.18% 

Models 



Length 

-1.000 

Thi'/Wt 

1.0533 

1.0555 

-0.21% 






GLOW 

0.9936 

0.9906 

0,30% 




Angle 

lOT) 

Thrust 

0.9965 

0.9956 

0.08%' 

Krigmg 

43 

4.67 

Height 

-0.873 

Weight 

0.7725 

0.7443 

3.79% 

ModeAs 



Length 

-1.000 

Thi'/Wt 

1.0506 

1.0568 

-0.59% 






GLOW 

0.9824 

0.9914 

-0.90% 

Minimize GLOW 




Angle 

0.616 

Thrust 

1.0013 

0.9957 


RS 

30.67 

3.33 

Height 

-1.000 

Weight 

0.8969 

0.8617 

4,09% 

ModfAs 



Length 

1.000 

Tlu'/Wt 

1.0251 

1.0286 

-0,34% 






GLOW 

0.9660 

1.0146 

-4.79% 




Angle 

(j:764“ 

Thrust 

1.0009 

1.0006 

0.04% 

Kriging 

57.67 

6.33 

Height 

-0.833 

Weight 

0.9060 

0.8732 

3.75% 

ModeAs 



Length 

0.676 

Thr/Wt 

1.0228 

1.0302 

-0.72% 






GLOW 

0.9675 

0.9680 

-0.05% 


Maximize Thriist /Weight Ratio 





Angle 


Thrust 

i.Ml6 


(1.57% 

RS 

27 

1 4 

Height 

-0.433 

Weight 

0.9450 

0.9073 

4.16% 
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Modfda 



Length 

1.000 

Thr/Wt 

GLOW 

1.0141 

0.9724 

1.0173 

1.0228 i 

-0.31% 

-4.93% 

Krig-ing 

Modeh 

62 

5 

Angle 

Height 

Length 

0.656 

-0.627 

1.000 

Thi'ust 

Weight 

Thr/Wt 

GLOW 

l.OOlC 

0.9385 

1.0157 

0.9690 

1.0014 

0.906:1 

1.0231 

0.9666 

om 

3.56% 

-0.73% 

0.25% 


*i4( + ) err<yrtermindicatesthaUhf:7nodeMsovf:7y7'eAlicting,a{-)indicaiesthatitisunde7'pi'edicting. 


The following obsen^at ions can Ix^ made based on the data in TABLE 5. 

• Avcragemirnber of arialysis arid gradient calls: In generah the RS models require fewer analysis and 
gradient calls to achieve the optimum than the kriging models do. This can Ix' attributed to the fact 
that the RS models are simple second order polynomials whereas th(' kriging iiKKlels are more 
complex. 

• Convergence rates: Although not shown in the table, optimization using the RS models tends to 

converge more quickly than wlien using kriging models. This can be inferred from the average 
number of gradient calls which is one to three calls fewer for the RS models. 

• Optiiminidesigiis: The optimum designs ol^tained from the RS and kriging nxxlels are essentially the 
same for each objective function, indicating that both ap})roximations send tlu^ optimization 
algorithm in the same general direction. The largest, discrepancy is the length for the minimize 
CLOW optimization; RS models predict the optimum GLOW occurs at the upjx^r lx>und on length 
(+1) while the kriging models yield 0.676. This difference is evident in FiC. 11. 

• To check the accuracy of tlu^ predicted optima, the optimum 
design values for angle, height, and length are used as inputs into the original analysis c<xles and the 
jXTcentag(^ difference between the ac'tiial and predicted values is computed. Tlu' pnxiiOion error is 
less than 5% for all cases and is 0.5% or less in three quarters of the results. 

In summary, then, neither model coiLsistently outperforms the other, and the difference in predictiv<^ 
capability of each model for each objective function is quite small. JliekTiyingiriodelsperforinaswellasnie 
sc cmidordc r RS mode Is e iw nt hough t h e gl obal portion o ft he krigingmode I i s only a cmis t an i . O ngoi ng work to 
improve incxlel accuracy and checking adequacy of fit is detailed in the next sect ion. 

6. Closing Remarks and Future Work. This work represents a preliminary investigation into 
the use of kriging as an alternative statistical-based approximation technique for modeling non-random, 
deterministic' computer experiments. A three variable engineering design example is used to compare the 
approximation capability of response surface modeling and kriging. The example is the multidisciplinary 
design of an aerospike nozzle which includes a CFD and a finite element model. With this simple, yet 
realistic engineering example, the use of kriging models as an alternative approximation technique has 
been demonstrated. At this point, there is inconclusive evidence to state that one aj)proximation method 
is more advantageous than another; however, the kriging models, using only a constant underlying global 
model and a Gaussian correlation function, perform as well as the sec'ond order response surface models. 

There are several research issues to address for the application of kriging and DACE metluxls for 
other (and larger) engineering design problems. 

• Fittingakrigirigjnodel: Fitting a kriging model requires solution of an k-dimensional, unconstrained, 
non-linear optimization problem, Eqn. (12), in order to determine the maximum likelihood estimates 
of the 0 parameters for the ''best'" kriging model. Pattern search methods and simulated annealing 
algorithms are currently being employed to jx^rform this optimization. For small problems with 
relatively few sample points, this optimization is rather trivial. However, as the size of the problem 
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increases and the number of sample points increases, the added effort needed to obtain the ^‘best'" 
kriging model may quickly begin to outweigh the iDenefit of building the approximation. 

• In this example, a constant is used for the global jwrtion of the kriging 
ino<lel based on the success of the work in [30j, [37j, and [47]. However, using a global polynomial 
model for f(x) in Eqii, (5) may further improve the accuracy of the kriging model. Giunta llj 
jx^rforms a preliminary investigation of such an approach and finds that minimal improvement was 
obtained. 

• Unlike RS model prediction, prediction with a kriging model requires 
the inversion and multiplication of several matrices: these matrices grow as the number of sample 
jx>ints increases. For large problems, prediction with the kriging model may l)ecome computationally 
exjx'iisive in and of itself. Furthermore, it is more difficult to look at a kriging model and determine 
the effects of the design variables on the response(s) since the global model is usually taken as a 
constant, and each prediction jx^int is essentially the sum of exponentially decaying functions based 
on R. 

• Validatiiigakrigh^^^ With RS models, R^ values and residual plots can Ix^ used to assess model 

fit and accuracy. Since kriging models interpolate the data, there are no residuals and alternative 
c hecks must lx‘ implemented to validate the mcxlel. In this example an additional 25 random data 
jx)ints are used to check mcxlel adequacy; however, more formal approaches exist. Otto, et al. [31] 
and ^32i and Yesilyurt and Patera [48] have develoix^d a Bayesian- validated surrogate approach which 
use's additional validation points to make qualitative assessments of the quality of the approximation 
mcxlel and provide theoretical lx)unds on the largest discrepancy between the mcxlel and the actual 
computer analysis. An alternative methcxl which clexjs not require additional ix)ints is Ic^ave-one-out 
cross validation [26]. Each sample point used to fit the mcxlel is removed one at a time, the model is 
rebuilt without the sample jx)int, and the difference l>etween the mcxlel without the sample jx>int and 
actual value at the sample point is computed for all of the sample points. Neither approach was 
implemented in this example due to the increased computation effort required. 

• Dtsupi of cxperbneiits for buildiiig kTigmg models: Are there designs which are better suited for 

sampling computer exix'riments and building kriging mcxiels than for sampling physical expe^riments 
and building RS mcxiels? The opinions on the appropriate experimental design vary; the only 
consensus reached thus far is that designs for non-random, deterministic computer analyses should be 
space filling. In this example, orthogonal arrays are used to build approximation following the work 
by Booker, et al, [3] and [4]. Giunta [11] uses D-optimal designs to fit hLs kriging and RS mcxiels; 
Sacks, et al. [38] suggest using IMSE-optiinal designs; and Koehler and Owen [20] discuss minimax, 
maximin, Latin hypercube, and scrambled net desigiLs for computer experiments. 

Despite the added complexity of fitting, using, and validating a kriging mcxlel, the px>tential gains in 
mcxlel accuracy justify continued investigation into the approach. The kriging software under 
development will fac-ilitate the use and validation of kriging mcxiels, increasing their attractiveness for 
enginc'ering applications. Finally, future work on the aerospike nozzle design problem includes expanding 
the scope' of the problem to include more design variables and responses and investigating the impact of 
decom|X)sing the problem into its disciplines by building approximation mcxiels of each discipline 
separately and examining the effect of different multidisciplinary design formulations (e.g., [l] and [8]) on 
the solution. 
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